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ABSTRACT 

We study the long-term dynamics of the PSR 1257+12 planetary system. Using the recently determined 
accurate initial condition by Konacki & Wolszczan (2003) who derived the orbital inclinations and the absolute 
masses of the planets B and C, we investigate the system stability by long-term, 1 Gyr direct integrations. 
No secular changes of the semi-major axes, eccentricities and inclinations appear during such an interval. 
This stable behavior is confirmed with the fast indicator MEGNO. The analysis of the orbital stability in the 
neighborhood of the nominal initial condition reveals that the PSR 1257+12 system is localized in a wide stable 
region of the phase space but close to a few weak 2 and 3-body mean motion resonances. The long term stability 
is additionally confirmed by a negligible exchange of the Angular Momentum Deficit between the innermost 
planet A and the pair of the outer planets B and C. An important feature of the system that helps sustain the 
stability is the secular apsidal resonance (S AR) between the planets B and C with the center of libration about 
180°. We also find useful limits on the elements of the innermost planet A which are otherwise unconstrained 
by the observations. Specifically, we find that the line of nodes of the planet A cannot be separated by more than 
about ±60° from the nodes of the bigger companions B and C. This limits the relative inclination of the orbit of 
the planet A to the mean orbital plane of the planets B and C to moderate values. We also perform a preliminary 
study of the short-term dynamics of massless particles in the system. We find that a relatively extended stable 
zone exists between the planets A and B. Beyond the planet C, the stable zone appears already at distances 
0.5 AU from the parent star. For moderately low eccentricities, beyond 1 AU, the motion of massless particles 
does not suffer from strong instabilities and this zone is basically stable, independent on the inclinations of the 
orbits of the test particles to the mean orbital plane of the system. It is an encouraging result supporting the 
search for a putative dust disk or a Kuiper belt, especially with the SIRTF mission. 

Subject headings: celestial mechanics, stellar dynamics — methods: numerical, N-body simulations — planetary 
systems — stars: individual (PSR 1257+12) 



1. INTRODUCTION 

Up to date, among about 130 extrasolar planetary sys- 
tems, the one around PSR 1257+ 12 remains the only one 
which contains Earth- sized planets t Wjals^cz^ n & Frail 19921 
IWolszczanl 1 1 994t lKonacki&~ W olszczanl l2003l) . It has been 
discovered with the pulsar timing technique that relies on ex- 
tremely precise measurements of the times of arrival (TOA) 
of pulsar pulses. Such a technique in principle allows a de- 
tection of companions as small as asteroids. In the case of 
the PSR B 1257+ 12 it enabled the detection of three planets 
A, B and C with the orbital periods of 25, 66 and 98 days and 
the masses in the Earth-size regime (see Table 1). Luckily, 
the two larger planets have the mean motions close to the 3:2 
commensurability which result in observable deviatio ns from 
a simple Keplerian descrip t ion of the motion (tRasio et al.l 
19921 iMalhotra et all ^991 iPeald Il993t IWolszczanl 1 19941 
Konacki et al. 1999). Konacki & Wolszczan] (120031) applied 
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the secular orbital theory of the PSR 1257+12 system from 
Konacki et al. (2000) to determine the masses and absolute 
inclinations of the orbits of the planets B and C. Recently, 
the same idea of incorporating the effects of mutual interac- 
tions between planets to remove the Doppler radial veloc- 
ity (RV) signal degeneracy present in the /V-Keplerian or- 
bital model (i.e., the undeterminancy of the system inclina- 
tion and the relative inclination of the orbits) has also been ap- 
plied to stud y the stro ngly interactin g, resonant system around 
Glies e 876 tLau ghlin & Chambers 200 jj iRivera & LissaueJ 
120011) . However, the accuracy of the timing observations by 
far exceeds the precision of the RV measurements. In con- 
seque nce, the accuracy of th e PSR 1257+12 fit obtained by 
Konacki & Wolszczan (2003) and hence the initial condition 
is superior to any initial condition derived from the fits to RV 
observations of solar-type stars with planets. Thus, the sys- 
tem constitutes a particularly convenient and interesting sub- 
ject for the studies of its dynamics. 

There have been few and limited attempts to determine the 
long-term stability of the PSR 1257+12 system. U sing the res- 
onance overlap criterion and direct integrations, iRasio et 
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(QUI, Malho tra et all 11991 and iMalhotral (Q993b) have 
found that the system disrupts on 10 5 yr timescale if 
the masses exceed abou t 2-3 masses of Jupiter. Also 
iMalhotra et alJ (119921) and Malhotra ( 1993a) have estimated 
that if the masses were about 20 — 40 M ffi then the system 
would be locked in the exact 3:2 resonances which would 
lead to the TOA signal very different from what is actually 
observed. In this paper, we investigate the stability of the 
PSR 1257+12 system in the Gyr-scale. We also perform a 
dynamical comparison of the PSR 1257+12 system to the in- 
ner Solar System (ISS) as the dynamics of the ISS is mostly 
driven by interactions between the telluric planets and is basi- 
cally decoupled from the dynamics of the outer Solar System 
( Lask arU994Hl997h . Finally, we carry out a preliminary anal- 
ysis of the stability of the orbits of massless particles (e.g. dust 
particles or Kuiper belt type objects) to determine the zones 
where such particles could survive and possibly be detected. 

2. NUMERICAL SETUP 

The stability analysis performed in this paper has mostly 
numerical character and is understood in terms of the maxi- 
mal Lyapunov Characteristic Number (LCN). Hence we treat 
the quasi-periodic orbital motions (LCN ~ 0) as stable and 
the chaotic ones (LCN > 0) as unstable. In order to re- 
solve the character of the orbits efficiently, we use the so- 
called fast indicator MEGNO (Mean Exponential Growth fac- 
tor of Nearby Orbits, Cincot ta & Sim6l l2000). The MEGNO 
is directly related to the LCN through the linear relation 
(Y) = at + b. For quasi-periodic motions a ~ and b ~ 2 
while for chaotic solutions a ~ (A./2) and b ~ where A, 
is the LCN of the orbit. However, there is no general rela- 
tion between the magnitude of the LCN and the macroscopic 
chang es of the orbital elements ( e.g., Murison et al.l 119941 
Micht chenko & Ferraz-Mellol 1200 111 . Even if the LCN is 
large, the system can stay bounded for a very long time. More 
insight into the relation of the degree of the irregularity of 
motion and the macroscopic ch anges of the orbits one can ob- 
tain with the FFT techniques ( Michtchenko & Ferraz-Mello 
120011) and its refined variant , the M odified Fourier Transform 
(MFT) developed bv lLaskarl (119931) . The MFT makes it possi- 
ble to resolve the fundamental frequencies in a planetary sys- 
tem a nd to determine their diffusion rates ( Robutel & Laskar 
1200 11) . In our tests, we use the MFT to verify the MEGNO in- 
tegrations for massless particles (Section 4) and to identify the 
mean motion resonances (MMRs). During the computations 
of MEGNO, for every planet k =A,B,C, and also for massless 
particles (k — 0) the complex functions fk(ti) = flA-expiA^ (?,-), 
are computed at discrete times f,-, where f, + i — ~ 10 days, 
over the time 2T between 13,000 P c and 64,000 P c (where 
Pc is the orbital period of the outermost planet C). Here, a\ 
denotes the semi-major axis of the relevant planet and A,£ is its 
mean longitude. In general, for a quasi-periodic solution of a 
planetary system, the frequency corresponding to the largest 
amplitude, a\, is one of the fundamental f requencies of mo- 
tion, called the proper mean motion, (Robut eT& Laskar] 
1200 11) . The 2-body MMRs can be identified through the con- 
dition qVj — pVj ~ 0, / 7^ j = A,B,C,0, where p,q > are 
prime integers. The MFT code gives us the rates Vi/Vj ~ p/q 
and then p and q are found by the continued fraction algo- 
rithm. Such resonances will be labeled as Pig : ?2p where Pi 
and P2 denote the bodies involved in the resonance. In this 
work we use the code which merges the MFT and MEGNO 
methods. It incorporates the MFT code which has been kindly 



provided by David Nesvorny 4 . 

If a planetary system is strongly interacting, MEGNO en- 
ables us to detect chaotic behavior typically over only 10 3 - 
10 4 orbital periods of the outermost planet. This feature is vi- 
tal for examining large sets of the initial conditions (see e.g ., 
Gozd ziewski et all200 It [Gozdziewski & Macieiewski 20031) . 
However, due to very small masses in the PSR 1257+12 
system, the mutual perturbations between the planets are 
very week. In order to resolve the character of such 
orbits, much longer integration times are required. We 
have determined that the general purpose integrators (like 
the Bulirsh-Stoer-Gragg, BSG) used to compute MEGNO 
are not efficient enough. Therefore, we have developed 
a sympl ectic scheme of com puting MEGNO ba sed on the 
idea o f Mikkola"& Innanenl Jl999) (see also iGozdziewskil 
(2003) and ICincotta et all 12*003)). Instead of directly 
solving the variational equations of the perturbed Kepler 
problem, one can differenti ate t he symplectic leap-frog 
scheme ( Wisdom & Holman 1991) and compute the varia- 
tions (and MEGNO) using the obtained symplectic tangent 
map. To use such a technique, a system of canonical vari- 
ables is required. In these variables, the planetary problem 
can be divided into integrable Keplerian motions and a per- 
turbation which is integrable in the absence of the Keple- 
rian part (usually it depends only on the coordinates). Our 
MEG NO code works in ternally in the Poincare coordinates 
(Laskar & Robutel 1995). As the integrato r core, the second- 
(LR2) and third- (LR3) order schemes by Laska rlfc Robutell 
(2001) are used. These integrators are much faster than the 
BSG algorithm. We can efficiently perform intensive com- 
putations of MEGNO (~ 10 4 points per stability map) for 
~ 10 6 Pc, keeping the fractional error of the total energy and 
the angular momentum below 10~ 10 when the 4 d step (for 
the LR3 integrator) is used. The tangent map approach is 
particularly suitable for the pulsar system due to typically 
small eccentricities of the investigated configurations. The 
method will be described in detail elsewhere (Breiter and 
Gozdziewski). 

However, one should be aware that even with such an ex- 
tended integration time-scale, in general we can only detect 
the strongest chaotic behavior originating mostly from the 
mean motion resonances between the planets. Investigations 
of subtle chaos originating in the secular system would re- 
quire much longer integration times, of the order of thousands 
of secular periods. Because these periods are in the range 
10 3 -10 5 yr (see Sect. 3) or longer, the total integration time 
should be counted in tens of Myr. In such a case, due to the 
computational limitations, direct approach applied in this pa- 
per is no longer practically possible. 

In the numerical experiments, we use the initial condition 
given in Tabled I t represents one of the tw o orbital configura- 
tions obtained by Konacki & Wolszczan (2003) transformed 
to the classical astrocentric elements. In some of the experi- 
ments these elements have been transformed to the canonical 
heliocentric elements inferred from the Poincare coordinates. 
In the second orbital configuration, the inclinations are such 
that /b,c — * 180° — /b,c- We assume that z'a obeys the same 
rule in that solution. It follows that both these solutions are 
dynamically equivalent and this reflects the fact that using the 
secular theory as the model of TOA measurements we are still 
not able to determine the absolute direction of the angular 
momentum vector. The TOA residuals exclude opposite di- 

4 www .boulder . swri . edu/ "davidn/fmf t/ 
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rections of the outermost orbits ( K onacki & Wolszczanl2003l) 
hence the orbits are either prograde or retrograde. 

3. STABILITY OF THE PSR 1257+12 PLANETARY SYSTEM 

We computed the MEGNO signature of the PSR 1257+12 
system over a few Myr (5 • 10 6 Pq). This test was repeated for 
many different, randomly selected initial variational vectors. 
The results of a few of the runs are shown in Figs^e. Both 
the regular evolution of MEGNO, Y(t), as well as the quick 
convergence of its mean value, (Y)(t), to 2, indicate that the 
system is very close to a quasi-periodic motion. Let us note, 
that our preliminary computations of MEGNO with the BSG 
integrator revealed a very slow divergence of MEGNO after 

0. 5 • 10 6 Pc from the theoretical value of 2. It can be explained 
by an accumulation of the numerical errors because this effect 
is absent in much longer runs employing the symplectic code. 
Even if such a weak instability were to be understood in terms 
of the chaotic behavior, it would correspond to an extremely 
long Lyapunov time, Ti = l/X, comparable to 10 9 yr (esti- 
mated through the fit of the relation (Y) = (X/2)t + b over 
t ~3.5-10 6 yr). 

In Fig. [2 we present the time-evolution of the canonical 
heliocentric elements. Clearly, the motion of the planets B 
and C is tightly coupled. This is illustrated in panels b and 
c which show the variations of the eccentricities and inclina- 
tions. Fig.^1 demonstrates the presence of a deep, so-called 
secular apsidal resonance (SAR), described by the argument 
= 03b — 03c (where 03 = co + Q. is the longitude of perias- 
tron; Q. and CO are the longitude of ascending node and the ar- 
gument of periastron of the planet, respectively), with a small 
semi-amplitude of the librations, ~ 50°, an d the apsides on av- 
erage a nti-aligned. In the discovery paper, Wolszczan & Frail 
( 1992) noticed the anti- alignment of the lines of apsides while 
iRasio et all (119921) and Malho tra etai] (f.1992) performed first 
theoretical explorations of the resonance in the framework of 
the Laplace-Lagrange theory. 

The SAR was recently found in many extrasolar systems 
discovered by the RV measurements. It is widely believed 
that the SAR, as the orbital state of multi-planetary sys- 
tems involving Jupiter-size planets, is crucial for maintain- 
ing th eir long term stability (for an overview see, LTi et alJ 
2003). However in certain cases, the SAR should be un- 
derstood as a typical feature of a planetary system when the 
secular angle 9 oscillates about 0° or 180°. In particular, 
this concerns the pulsar system: an inspection of Table 3 
shows that the proper frequencies of the pericenter motion, 
g p , are not related through any simp le linear comb ination ful- 
filling th e usual resonanc e r ules. iLee & Peald C2003) and 
Micht chenko & Mamotral (120041) show that the orbital config- 
urations of two planets are generically stable if the system is 
far from strong mean motion resonances and collision zones, 

1. e., when the assumptions of the averaging theorem are ful- 
filled. In such a case, the secular character of the system de- 
pends on the variations of the angle 9. Three different modes 
of its motion are possible: circulation, oscillations around 0° 
or 180° and oscillations about 0° in the regime of large eccen- 
tricities. The first two modes are know n from the classical 
Lapla ce-Lagrange linear secular theory (Murrav & Dermott 
2000) while the third one was discovered by the cited authors. 
The third mode corresponds to the true, non-linear resonance 
and the secular system may be unstable in its neighborhood. It 
should be noted that for the two first modes there is no zero- 
frequency separatrix between the circulation and oscillation 
regimes. In this sense, such a SAR is not a true resonance. 



Due to a very small mass of the planet A, the pulsar system 
may be approximated by a 3-body model with equal planetary 
masses, s mall e ccentricities and inclinations. In such a case, 
iPauwelsl ([E>83) showed that, at least in the linear approxi- 
mation, the SAR regime covers almost the entire phase space 
of the secular planetary system, the occurrence of a SAR is 
almost a certainty and the system is stable. This could ex- 
plain the stability of the pulsar system. However, some of 
its specific dynamical features break the assumptions of the 
secular model. Namely, the proximity of the system to the 
first-order 3:2 MMR and other mean motion resonances. In 
order to estimate the influence of the near 3:2 MMR on the 
secular solution, we integrated the system over 10 Myr sam- 
pling the orbital elements every 100 yr. Subsequently, the 
MFT was applied to analyze the complex signals e p exp(i03 p ) 
and sin(/p/2)exp(i£2 p ) derived from the heliocentric canoni- 
cal elements. In this way, we can obtain an estimate of the 
precessional frequencies of the periastrons, g p , and the nodal 
lines, s p . These frequencies, together with the proper mean 
motions, n p , compose the set of fundamental frequencies of 
the pulsar system. Their values and the corresponding periods 
are given in Table 3. The periastron frequencies, measured in 
arcsec/yr, are ~ 197.742, ~ 43.881 and ~ 13.157. The first 
one is substantially different from its analytical approxima- 
tion (in terms of the Laplace-Lagrange theory). This will lead 
to a fast discrepancy between the analytical and numerical so- 
lutions. Clearly, the classical secular theory has to be modi- 
fied in order to account for the effects of the near 3:2 M MR. 
Such a theory has been already developed by Malhot ra et alJ 
(1989) for the satellites of Uranus. It can also be applied to 
the PSR 1257+12 system however it is somewhat beyond the 
scope of this paper. 

Usi ng the symplectic integrator WHM (Wisdom & Holman 
119911) from the SWIFT package llLevison & Duncan! 1994T as 
well as the LR2 integrator, we also performed a few long term, 
1 Gyr, integrations of the PSR 1257+12 system. In the first 
case, the time step equal to 1 d resulted in the fractional error 
of the integrals of the angular momentum and the total energy 
at the level of ~ 10~ 10 . The LR2 algorithm with the 4-day 
integration step resulted in even higher accuracy. During the 
1 Gyr interval, no signs of instability are observed. We have 
not detected any secular changes of the semi-major axes, ec- 
centricities and inclinations. Their values stay within the lim- 
its shown in Fig.^and the SARpersists with an unchanged 
amplitude of the librations (Fig. |2ji) over 1 Gyr. 

Lest us note that from the 1 Gyr integrations it follows that 
the orbital elements of of the innermost planet A vary in a 
regular way (Fig. We can also observe a time evolution 
of the argument 9i = 03a — 03b as a semi-regular sequence 
of rotations alternating with irregular "librations" about 180°. 
This effect is preserved over the entire period of 1 Gyr (see 
Fig. |2j;,d) and may indicate that the orbit crosses the sepa- 
ratrix of a resonance. This is illustrated in Figs |2j;,d which 
show the eccentricity of the planet A (multiplied by 10,000), 
<?a, plotted together with the argument 9i. Clearly, there is a 
correlation between the librations of 9i and small values of 
ca- We also show ec(9) and eA(9i) collected over the 1 Gyr 
integration (Fig. Et,f). The first case represents a trajectory 
in the resonant island of the SAR whereas there is no clear 
sign of librations in the eA(9i) plot. In fact, this effect is 
only geometrical in nature and can be explained by estimat- 
ing the secular frequencies of the system using the well known 
Lagrange-Laplace theory. If the MMRs are absent and the dis- 
turbing function is expanded up to the first order in the masses 



4 



Gozdziewski, Konacki & Wolszczan 



and to the second order in the eccentricities and inclinations, 
the equations of the secular motion are integrable. Their so- 
lution, relative to the eccentricities and the longitudes of peri- 
astron, is given in terms of the so called eccentricity vectors, 
(h,k) = (e sin 03, e cos 03), by (Murr ay & Dermotl2 000) 

hj{t)=Y t ej l iwa(g i t + ^ i ) 

i 

M0=L e -'V cos (s j *+P'')> 

i 

for every planet j =A,B,C. The constant amplitudes e ; - ,■ and 
secular frequencies gi are the eigenvectors and eigenvalues 
of a matrix with coefficients given explicitly in terms of the 
masses and constant semi-major axes of the planets. The scal- 
ing factors for the eigenvectors <?, ; and the phases p, are deter- 
mined by the initial condition. Geometrically, the time evo- 
lution of every eccentricity vector can be describe d as a su- 
perpos ition of the eigenmodes corresponding to gi l lMalhotral 
1993b). The parameters obtained for the PSR 1257+12 sys- 
te m are given in Tab le El They are consiste nt with the results 
of lMalhotral(ll993al) and lRasio et alJ i 19921) who analyzed the 
PSR 1257+12 system involving the two larger companions. 
They found that the secular evolution of the eccentricity vec- 
tors of the two outer planets is almost entirely driven by the 
eigenstate corresponding to g\. Since we have eg a — —ec\ 
and the other components of ej are much smaller and almost 
equal, it follows that (Ob — COc + 180°, which corresponds to 
the SAR of the planets B and C. For the planet A, the eccen- 
tricity vector is a superposition of all the three eigenmodes 
with the leading second and third mode having comparable 
amplitudes. It appears that using the secular approximation, 
we can explain the semi-librations of 8i, The minima of 
occur when the modes corresponding to the frequencies g2 
and gi are in anti-phase. Because the amplitudes eA.i and eA,3 
have the same sign and similar magnitudes while e^.i has a 
much smaller magnitude, the following condition has to be 
satisfied to grant a minimum of e& 

(S2f + P 2 )- (^ + (33)^180°. 

Since p 2 — P3 ~ -180° (see TableE), we have {g%- gi)t = 
360° and the period of anti-alignment is P23 = 360°/(g2 — 
#3) ~ 43,000yr. Curiously, from the relation between gi, g\ — 
7(g2 — gj) — 0.001°/yr, it follows that this period is almost 
commensurate with Pi (where Pi = 360° /g\). 

Near the moment of anti-alignment, the eccentricity of the 
planet A, e^, is driven mostly by the gi-mode that is in anti- 
phase with the eccentricity vector for eg. This leads to the 
quasi-librations seen in Fig. E] as they repeat with the same 
period as the variations of e&, P23. 

3.1. Comparison to the inner Solar System 

The results of the 1 Gyr integration, expressed in the canon- 
ical heliocentric elements, enable us to analyze the time- 
evolution of t he so-called Angular Momentum Deficit (AMD, 
lLaska.11 9971): 

C = L ,r + m VAV^ 1 - V 1 - g | cos 'p). 

p=A,B,C m P + M * 

where p. p = G(m p + m*), G is the gravitational constant, 
a p ,e p ,i p is the semi-major axis, eccentricity and inclination of 
a planet relative to the invariant plane and is the mass of the 
central body. The AMD indicates the deviation of planetary 
orbits from a stable circular and coplanar motion for which 



C = 0. This quantity is preserved by the averaged equations of 
motion (Laskar 1997) and its stability provides the stability of 
the secular system in the absence of short-period resonances. 
The AMD can be understood as the amplitude of the irregu- 
larity present in the averaged system. Large values of AMD 
lead to a chaotic motion and for a certain critical value to a 
crossing of th e orbit s and the disruption of a planetary s ystem 
(Laskar 1997, 2000; Michtch enko & Ferraz-Mellol2001l) . 

The PSR 1257+12 system is close to the 3:2 commensura- 
bility. Its critical argument circulates but it does not mean 
that the time-averaged effects of the n ear-resonance vanish 
(Malh otra et al.lll989tlMalhotralll993al) . Hence, the applica- 
bility of the AMD signature in the studies of the stability in the 
real, unaveraged system can be problematic. Instead, the inte- 
gral obtaine d by aver aging th e quasi-resonant sy stem should 
be applied (Michtchenk o & Ferraz-Mellol 120011) . Neverthe- 
less, we have decided to calculate the AMD and to examine 
its behavior in order to compare the results with those ob- 
tained for the Solar sy stem (which contains t wo planets close 
to 5:2 MMR, see e.g.. llto & Tanikawall2002l) . The time evo- 
lution of AMD in the pulsar system is shown in Fig. [2] The 
AMD stays well bounded and very regular. There is very lit- 
tle exchange of the AMD between the innermost planet A and 
the pair of the bigger companions B and C. It suggests that 
the motion of the planet A is decoupled from the dynamics of 
the planets B and C in the long-term scale. This situation is 
qualitatively different from the inner Solar system (ISS). The 
AMD of the ISS is not strictly preserved due to the pertur- 
bations of Jupiter and Saturn, nevertheless it can be consid - 
ered roughly constant tLaskaJI 19971 llto & Tanikawal 120021) . 
Ito & Tanikawa (2002) published the results of 5 Gr integra- 
tions of the ISS which reveal rapid AMD variations of Mer- 
cury. They are much larger than the changes of AMD for the 
Venus-Earth pair. The variations of the AMD of Venus, Earth 
and Mars are also irregular and substantial. In fact, t he ISS 
is cha o tic having the Lyapunov time of about 5 Myr (Laskar 
1994). lLaskarl (11994) found that this chaos is physically sig- 
nificant as it can lead to the ejection of Mercury from the in- 
ner ISS during a few Gyr. Howev er, the source of c haos in 
the ISS is still not well understood dLecar et all2001l) . In this 
sense, the character of the motion of the PSR 1257+12 system 
is quite different. The evolution of its AMD is very regular 
in spite of smaller distances, larger masses and thus stronger 
mutual interactions between the planets. Unlike Mercury's, 
the AMD of the planet A is negligible when compared to the 
AMD of the B-C pair as it contributes only about 1/1000 of 
the total value. On the other hand, the orbital coup ling of the 
pairs B-C and Venus-Earth dlto & Tanikawa 2003) is a sim- 
ilar feature of the inner Solar and PSR 1257+12 systems al- 
thoug h it has a d i fferen t nature. The long term integrations by 
llto & T anikawa (2002) revealed an anti-correlation between 
the changes of the orbital energies of Venus and Earth and, 
simultaneously, a correlation between the changes of their 
orbital angular momenta and the eccentricities. These ef- 
fects can be explained through the influence of the Jovian 
planets. Hence their dynamical source is external for the 
coupled planets while the coupling of the B-C pair in the 
PSR 1257+12 system is provided by the anti-aligned SAR. 
Finally, the chaotic orbital evolution of the ISS may signif- 
icantly depend on the weak coupling with the outer planets. 
An equivalent effect is obviously absent in the PSR 1257+12 
system. 
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3.2. Dynamical environment of the PSR 1257+12 planetary 

system 

In the next set of experiments, we look at the initial condi- 
tion in a global manner in order to find out whether the cur- 
rent state of the system is robust to the changes of the formal 
initial condition (IC). However, thanks to a very precise deter- 
mination of the initial condition such changes, if considered 
consistent with the TOA measurements, are very limited. For 
example, the formal l g error of the semi-major axes inferred 
from the parameter x° (Konacki & Wolszczan 2003) is at the 
level of 10~ 6 AU! Nevertheless, the localization of the IC in 
the phase space (e.g. its proximity to unstable regions) is crit- 
ical to verify its character. The examination of one isolated IC 
does not provide a definitive answer to the question of stabil- 
ity. 

3.2.1. Mean motion resonances 

We computed one-dimensional scans of MEGNO, (Y), by 
changing the semi-major axis of one planet and keeping 
the other orbital parameters fixed at the values given in Ta- 
bleQ The results of this experiment are illustrated in Fig. [3] 
The MEGNO scans were computed with the resolution of 
~ 5 • 10~ 6 AU and 10~ 6 AU for close-up scans. They re- 
veal a large number of spikes, some of them very close to the 
nominal positions of the planets marked with the large filled 
dots. Most of these spikes can be identified as 2-body MMRs 
between the planets. Currently, the IC is well separated from 
low-order MMRs. But some unstable high-order MMRs ap- 
pear close to the nominal IC. The most relevant ones seem to 
be C31:B21 MMR (i.e. between the planets B and C, see the 
bottom panel in Fig.|3j and C53:B36 MMR (shown in the scan 
for the planet B). Also the planet A is close to B29:A1 1 MMR. 
Yet, such resonances are extremely narrow. Their widths can 
be roughly estimated as less than 5 • 10~ 6 AU what follows 
from their shapes and the resolution of the MEGNO scans. It 
is unlikely that they can affect the regular motion of the sys- 
tem at its nominal position. However, altering the semi-major 
axis of the planets by small shifts, ~ 5 • 10~ 5 AU, would push 
the system into these unstable regimes. In order to demon- 
strate it, we calculated MEGNO for a configuration corre- 
sponding to C31:B21 MMR (see the right column of Fig.0). 
This resonance i s the inclinat ion-type MMR (according to the 
terminology by Peala [19761) . Clearly, its critical argument 
O = 3lXc~ 21A.B — 10£2b exhibits a sequence of alternating 
librations and circulation giving rise to the unstable behavior. 
Let us note, that the chaos is formally extremely strong. The 
Lyapunov time, estimated by the linear fit to the MEGNO, is 
only about 3000 yr (see Fig.0^). 

The 2-body MMRs cannot explain all the peaks of MEGNO 
present in the scans. Some of them seem to be related to the 3- 
body resonances involving all three planets. The importance 
of 3-b ody MMRs was studied by | N esvorny & Morbidelli 
( 119981): iMorbidelli & Nesvornvl ( 11999ft . iMurrav et alJ (119981) 
and [Murray & Holman (1999). They found that weak 3- 
body MMRs can strongly influence asteroidal motion and 
explain the short Lyapunov times which cannot be under- 
stood if only 2-body MM Rs of high-order are considered. 
IMurrav & Holmanl lll999) proved that the chaotic behavior of 
the Outer Solar system is governed by the overlapping 3-body 
MMRs involving Jupiter, Saturn, and Uranus. It is obviously 
interesting to verify if some of the instabilities visible in the 
a-scans can be related to such 3-body resonances. 

A 3-body resonance can be defined by the following condi- 



tion (iNesvornv & Morbidelh11 998): 

where A,p is the mean motion of the given planet. A criti- 
cal argument of such resonance is a linear combination of the 
longitudes, arguments of pericenters and nodal longitudes: 

a, A :, B :! C = 'A^-A + 'B^B + 'cA-C + 

PA® A + PB&B + PC&C + <7A^A + <?B^B + <7C^C 

whose integer coefficients i v ,p p ,q p fulfill the d'Alambert rule 
Lp=A,B,c('p + Pp + c lp) = and the usual requirement of rota- 
tional symmetry. Due to the small ratios of the secular fre- 
quencies to the mean motions in the PSR 1257+12 system, 
~ 10~ 5 , the resonance condition may be approximated by 
;'a«a + *b«b + 'c«c — 0. For example, one of the MEGNO 
peaks for ac — 0.46601 AU, can be explained as the com- 
bination of the proper mean motions 3«a — 14«b + 9«c — 
0.01°/yr. Hence it corresponds to the 3-body MMR 3:-14:9. 
It is separated from the the nominal position of the planet 
C by only ~ 3 • 10~ 5 AU which is at the 3o error level of 
ac- The time evolution of the critical argument, G3 : _i4 : 9, 
shown in the left column of Fig.|4l confirms the slow changes 
of this angle as well as a presence of quasi-librations alter- 
nating with circulations. This indicates a chaotic evolution. 
The MEGNO computed for the corresponding initial condi- 
tion does not grow fast, at least over a few Myr used in this 
test. This mean that the chaos is moderate. Nevertheless, the 
long term stability of the system (likely, the chaos would af- 
fect mostly the motion of the planet A) could be confirmed 
only by direct numerical integrations. Our 1 Gyr integrations 
for the ICs of the two MMRs near the planet C do not reveal 
any secular changes of the elements. Likely, in order to prop- 
erly investigate subtle effects of these resonances and their 
influence on the system, a refined model of the motion should 
be used. Such model should incorporate the effects of general 
relativity, a possible error in the pulsar mass and other factors, 
like the dependence on the observationally unconstrained or- 
bital elements of the planet A. 

Other low-order 3-body resonances that can be identified in 
the Fig.[5]are 1:-2:1, l:-4:2, -2:10:7 (almost overlapping with 
B71:A27), 3:-16:12, l:-6:5 (see the flA-scan). Finally, also an 
overlapping of 2-body MMRs is possible. One such example 
is marked in the ae-scan where in the neighborhood of the 
3:2 MMR, a sharp peak of MEGNO at the simultaneous po- 
sition of B60:A23 and C76:B51 MMRs is present. All these 
MMRs are very narrow since the width of 3-body MMRs is 
proportional to the masses in the second order and these 2- 
body MMRs are of high order. They are relatively distant 
from the nominal positions of the planets thus it is unlikely 
that they can affect the motion of the system. 

Finally, Fig.|4]shows the effect of varying the assumed mass 
of the neutron star on the determination of the semi-major axis 
of the outermost planet and the localization of the resonances. 
Obviously, the sequence of MMRs does not change. A dif- 
ferent mass of the pulsar can lead to a substantial shift of the 
initial semi-major axis (in our case, ac) and the weak unsta- 
ble resonances, like the C31:B21 MMR, may end up much 
closer to the nominal positions of the planets (for smaller than 
canonical mass of the parent star, M psl - = 1 .4M ). 

3.2.2. Stability in the (a p ,e p )- and (eg,ec)-planes 

The results of the 2-dimensional MEGNO analysis are 
shown in Fig. [6] In these maps we analyze the influence of 



6 



Gozdziewski, Konacki & Wolszczan 



the initial eccentricities of the planets on their motion. We 
used the initial conditions from Table 1 and changed those ele- 
ments that correspond to the coordinates in the stability maps. 
Due to a lower resolution (400 points in the semi-major axis 
range), some features present in the one-dimensional scans 
are not so clearly visible. Nevertheless, in the (g!b,<?b) -plane 
(the middle panel), the thin strips of the MMRs can be eas- 
ily identified. These maps additionally reveal the width of the 
strongest MMRs (like C3:B2, C13:B9 and C10:B7) and the 
border of the global instability. These stability maps confirm 
that the nominal positions of the two bigger planets, marked 
with filled circles, are far from strong instabilities of the mo- 
tion. Also in the (e%,ec) -plane, the nominal IC lies in a wide 
stable region, far from zones of chaotic motion (see Fig. 0. 
The (oa, <?A)map for the planet A shows a dense net of narrow 
unstable regions, some of which are close to the nominal po- 
sition of the planet. They have been already identified in the 
1 D scans for a a . Let us note that similar maps were computed 
bv lFerraz-Mello & MichtchenkoH2002l) for the planet B with 
their FFT fast indicator. Their results are generally consistent 
with ours even though they used a different initial condition 
of the pulsar system. 

During the MEGNO integrations we have also computed 
the maximum value of the SAR argument 9 = (53b — 05c with 
respect to the libration center of 180°. It enables us to estimate 
the semi-amplitude of librations and the extent of the SAR in 
the space of the scanned orbital parameters. The result of 
such an experiment conducted in the (ab , <?b ) plane are shown 
in Fig. [8] It uncovers an extended zone of the SAR with the 
smallest semi-amplitude of 9 in the vicinity of the nominal 
initial condition of the PSR 1257+12 system. The near 3:2 
MMR is also clearly present in this map. 

3.2.3. Limits on the unconstrained parameters of planet A 

So far, after Konacki & Wolszczan! (|2003), we were using 
the average inclination of the orbits of the planets B and C as 
the orbital inclination of the planet A. However, we can try 
to verify whether the dynamics can provide any limits on the 
inclination and hence mass of the innermost planet A. To this 
end, we computed MEGNO for the inclination of the planet 
A in the range [30°, 70°] and the position of its nodal line in 
the range [—90° : 90°], with the resolution of 1° in both coor- 
dinates. Thus we varied both the mass (to preserve mAsinz'A 
determined from the TOAs) as well as the relative inclina- 
tion of the planet A to the orbital planes of B and C. The 
result of this experiment is shown in Fig. |9] The MEGNO 
scan (the left panel of Fig. [9) reveals a well defined stable re- 
gion. The assumed position of the planet A is in the center 
of this zone. To obtain this picture, an extended integration 
time span of 10 6 fc was used. This choice was dictated by 
the analysis of a few orbits of A which are substantially in- 
clined to the orbital planes of B and C. An example corre- 
sponding to the initial £1a = — 75° (thus located in the unsta- 
ble region) is shown in Fig. [H)] For the integration time up 
to about 0.5 Myr (~ 5 ■ 10 5 P C ) MEGNO stays close to 2 but 
then it suddenly grows, indicating chaotic behavior. To ex- 
plain its source, we computed the eccentricity, inclination 
and the argument of pericenter, g&, of the planet A relative to 
the invariant plane. Evolution of these elements is shown in 
FigEll Clearly, both <?a and exhibit long-term, large am- 
plitude oscillations which are exactly in an anti-phase. The 
period of these oscillation is relatively very long, of the order 
of 10 4 yr. Simultaneously, g& temporarily librates about 90° 



or 270°. These librations indicate that the pericenter preces- 
sion of A stops. S uch features are typical for the well known 
Kozai resonance (lKozail ll962) found for highly inclined as- 
teroids in the Solar system. The argument of pericenter g\ 
can be treated as the critical argument of this resonance. Be- 
cause librations are followed by circulations of the critical ar- 
gument, the orbit crosses the separatrix which explains why 
the motion eventually becomes chaotic. 

The effects of the Kozai resonance are illustrated in the right 
panel of the e max -map (the maximal eccentricity attained by 
the planet A during the integration time) in Fig.|9] The sharply 
ending zone of moderate eccentricities is narrower than the 
stability region. We should be aware that in the transient ar- 
eas the integration time could still be to short to detect the 
irregular motion. Note that around the border of the stable 
zone, e\ can be as large as 0.6-0.8. This resonance puts sig- 
nificant limits on the position of the nodal line of the planet 
A as well as its relative inclination to the orbital planes of 
B and C. Assuming that the motion of all three planets is in 
the same direction, the nodal line of A cannot be separated 
by more than about ±60° from the nodes of B and C and the 
relative inclination should not exceed ~ 45°. Otherwise, the 
eccentricity of A would be excited to large values, > 0.7 over 
only 10 5 yr and for e\ ~ 0.8 the orbits of A and B would cross 
or close encounters would become possible. 

4. DYNAMICS OF MASSLESS PARTICLES 

While investigating the dynamics of extrasolar planetary 
systems one can ask whether the minor bodies can survive 
in the gravitational environment of the giant primary bodies. 
These can be small telluric planets in the systems with Jupiter- 
like planets. We can ask the same question about asteroids, 
cometary bodies or dust particles in the PSR 1257+12 sys- 
tem. Because the motion of the planets appears to be strictly 
stable, we can consider a simplified, restricted model. In this 
model one assumes that a probe mass moves in a gravita- 
tional tug of the primaries but does not influence their motion. 
In the simplest case such a model is we ll known as the re- 
stricte d, planar circular 3 -body problem ( M urray & Dermotfl 
2000). This simplification helps us to explore the phase space 
of the planetary system with reasonable computational effort. 
In this model, the MEGNO indicator and the MFT are evalu- 
ated only for the probe mass. In the numerical experiments, 
we investigate the orbital stability of massless particles in a 
few regions of the orbital space: between the planets A and B 
(region I), B and C (region II) and beyond the planet C (region 
III). 

In the first series of numerical tests, we vary the initial semi- 
major axis of the probe mass, ao, and its eccentricity, <?o, while 
the initial inclination is constant, z'o = 50°. Thus the tested 
orbits are slightly inclined to both orbital planes of the planets 
B (/ B = 53°) and C (z' c = 47°). The angular variables of the 
massless particles are set to Q.q — G3o = Mq = 0°. This way we 
essent ially follow the remarkable work of Robutel & Laskai 
(2001) who investigated the short-term dynamics of massless 
particles in the Solar System using the MFT. These authors 
argue that the global picture of the motion in the restricted 
problem does not qualitatively depend on the initial phases of 
the test particle. The initial plane (ao,eo) is representative for 
the dynamics because it crosses all resonances. By changing 
the initial orbital phases of the massless particles, the width 
of the resonances may vary but their influence on the motion 
can still be detected. 

In order to check the MEGNO signatures (e.g., to verify 
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whether the integration time is not too short), during the in- 
tegrations, we compute the diffusion rate, Go = 1 — rfi'/rft-', 
where m 1 ) and are the mean motion frequencies obtained 
over the intervals [0,T] and [T,2T], respectively, where T 
(which we will call the base period from hereafter) has been 
set accordingly to the investigated range of the semi-major 
axes. The close encounters with the planets are controlled. 
We assume that if the distance between the probe mass and 
a planet is less than the Hill radius, r# = [m p / (3m*)] ^/^a^, 
then the particle has collided with the planet or its motion be- 
come strongly chaotic. 

Within this numerical setup, in the region II we have not de- 
tected any stable motions. The region I turns out to be a much 
more interesting one. The results are shown in the MEGNO 
and o maps in Fig.^20 e ft panels for the region I, middle and 
right panels for the region III). The integrations have been 
carried out for the total time, 2T, of about 32,000 Pq- In both 
maps we can clearly recognize the planetary-crossing lines 
marked by large values of both indicators which correspond to 
strongly chaotic motions. These lines (white lines in Fig. II 1> 
are the solutions to the equations a p {\ +e p ) = ao(l — eo) 
for ao > a p and a p (\ —e p ) = ao(\ +<?o) for ao < a p where 
a p ,e p ,ao,eo denote the initial semi-major axis and the eccen- 
tricity of a planet and a test particle, respectively. 

In this area and above the collision lines, the particles are 
scattered mostly by close-encounters and collisions with the 
planets A and B. Remarkably, a stable area exists under the 
crossing-lines for the planets A and B. It is divided by a num- 
ber of MMRs with the planets. In order to verify this result, 
we have integrated the motion of 200 particles spread over 
[0.19,0.35] AU with the same resolution in the semi-major 
axis, ao, as in the 2D MEGNO test. These integrations were 
continued up to 15 Myr (here, we used the RMVS3 integra- 
tor from the SWIFT package). A resulting one dimensional 
scan of MEGNO over aq is shown in Fig. ^1 in the re- 
gions classified with MEGNO as regular, the massless par- 
ticles have mostly survived during the integration time while 
in the chaotic or close to chaotic areas, they have been quickly 
removed or collided with the planets (in this test the same cri- 
terion of 1 Hill radii for a collision event is used and in such 
a case we set logOo = 1). Some discrepancies are most likely 
related to the different integrators used in the experiment; the 
Bulirsh-Stoer integrator follows the orbits of massless par- 
ticles much more accurately that the symplectic integrator. 
This experiment enables us also to independently estimate the 
proper MEGNO integration time and to "calibrate" the scale 
of the diffusion rate, Oo- Comparing the MEGNO and Oq- 
maps, we conclude that for logo less than about (—7)— (—8) 
the motion can be considered as regular (also compare the 
maps in the left column of Fig. II 1> . Let us note that both indi- 
cators are in an excellent accord and the diffusion rate calcu- 
lated by the MFT algorithm seems to be even more sensitive 
to an unstable motion than the MEGNO is. 

Using a similar numerical setup, we have investigated the 
region III. However, due to an extended range of the semi- 
major axes, we have divided it into two parts: the region Ilia 
between 0.5 AU and 0.97 AU and the region Illb between 
0.9 AU and 3.9 AU (a distance comparable to the radius of the 
Asteroid Belt in the Solar System). The maps of MEGNO and 
the diffusion rate, logo, in the (ao,eo)-plane for the region 
Ilia (obtained after the time 2T ~ 32,000 Pc), are shown in 
the middle column of Fig. while the maps for the region 
Illb (2T ~ 64,000 Pc) are shown in the right column. For 



the region Ilia, in both maps, we can clearly recognize the 
collision line with the outermost planet while the crossing line 
with the planet B is in the zone of strong chaos. A number of 
MMRs appear as narrow vertical strips. This test shows that 
for a moderately low initial eccentricity, eo, the stable zone 
is extended and begins just beyond the orbit of the planet C. 
Obviously, with a growing eo the border of the stable zone 
shifts toward a larger ao (at the distances of about 1 AU the 
zone of stability reaches eo — 0.5). For the region Illb (ao S 
[0.9,3.9] AU), the results are shown in the right column of 
Fig. ^2 F° r efficiency reasons, in this test the step in eo is 0.1 
(the resolution in ao was left relatively high, at 0.005 AU). 
The collision lines are clearly present. Otherwise, this zone is 
mostly regular even for very large eo > 0.5. Similarly to the 
region Ilia, the scans reveal some extremely narrow unstable 
MMRs, most of them with the two outer planets. 

Finally, we examined the stability of orbits inclined to the 
mean orbital plane of the system. In this experiment we tested 
the motion of particles in the region EH. The results for the re- 
gion Ilia are shown in Fig.^](a number of additional scans, 
not shown here, allow us to extrapolate the results obtained for 
this region to the region Illb). In the first two experiments, we 
scanned the (ao,eo) space for the initial ;'o set to 75° and 87°, 
respectively. These inclinations, taken relative to the mean 
orbital plane of the PSR 1257+12 system, are comparable to 
the inclinations of the Kuiper belt objects in the Solar System. 
Again, the collision lines and the net of MMRs clearly appear 
in both the MEGNO and Oo scans which are shown in the the 
left and the middle panels of Fig. ^] The zone of strong chaos 
covers the crossing zones with the planets B and C while the 
collision line with the planet A is much more narrow and sep- 
arated by a quasi-regular area in w hich Op ~ 10~ 6 . Such an 
effect has been already observed bv lRobutel & Laskail J2001I) 
— for higher inclination the close encounters with the plan- 
ets are less frequent than for moderately small inclinations 
which explains the smaller extent of the strong chaos. More- 
over, after sufficiently long time the massless bodies will be 
removed from above the collision lines, excluding cases when 
the probes are trapped within stable MMRs. In the next test 
illustrated in the right panel of Fig.^J eo was set to and the 
initial inclination ;'o was varied in the wide range [10°, 90°]. 
This experiments allow us to generalize the results of the pre- 
vious scans. The stability of inclined particles is basically 
independent on the initial inclination, at least for moderately 
small eccentricities. 

One should be aware that the results of the above analysis 
should be treated as preliminary ones. Our study is restricted 
both to the short-term dynamics and a small part of the pos- 
sible volume of the parameter space. However, the timescale 
of the integrations is long enough to detect the most unstable 
regions in the phase space and to point out the regions where 
the particles can be long-term stable. Also in some chaotic 
regions the particles can still persist over a very long time but 
this can be verified only by direct integrations. Let us also 
note that in order to derive rigorous estimates of the proper 
mean motion frequencies one has to employ angle-action like 
coord inates, e.g., the Poincare coordinates ( Laska r^fe Robutell 
119951) or the Jacobi coordinates. Our integrations for mass- 
less particles were carried out in the Keplerian, astrocentric 
coordinates and the osculating elements analyzed by MFT are 
inferred from these non-canonical coordinates. Due to very 
small masses involved, the effects caused by the use of non- 
canonical coordinates are negligible. 
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5. CONCLUSIONS 

In this work we carry out numerical studies of the stability 
of the PS R 1257+12 system using the in itial condition deter- 
mined by Konacki & Wolszczan (2003). The long term in- 
tegrations utilizing the symplectic integrators, extended over 
1 Gyr, do not reveal any secular changes in the semi-major 
axes, eccentricities and inclinations of the planets. Using the 
notion of the Angular Momentum Deficit (AMD), we do not 
find any substantial exchange of the angular momentum be- 
tween the innermost planet and the pair of the outer, bigger 
planets B and C. The AMD of the planet A is negligible when 
compared to the AMD of the B-C pair. This is very different 
from the case of the inner Solar System in which the varia- 
tions of the AMD of Mercury are the most significant ones. 
The PSR 1257+12 system has the MEGNO signature typical 
for a strictly regular, quasi-periodic configuration. The two 
outer planets are close to the 3:2 mean motion resonance and 
are orbitally tightly coupled. The presence of the secular apsi- 
da l resonanc e is quite typical for suc h system as demon strated 
by Pauwels (1983) and recently by Lee & Peale (2003) and 
Micht chenko & Malhotral ( 12004) . The semi-amplitude of the 
critical argument is about 50° and it persists in a wide range of 
the orbital initial parameters. The SAR in the PSR 1257+12 
system i s yet an other such case among extrasolar planetary 
systems dJi et all2003l) . 

The neighborhood of the nominal initial condition is inves- 
tigated by calculating the MEGNO signature in a few rep- 
resentative regimes of the semi-major axes and eccentricities 
of the planets. These stability maps reveal that the nominal 
initial condition is located in an extended stable zone, rela- 
tively far from any strong instabilities of the motion. How- 
ever, numerous weak mean motion resonances can be found 
in close proximity to the nominal positions of the planets. 
These are both 2-body resonances between the planets (like 
the 31:21 MMR between B and C) and 3-body resonances, 
among which the 3: -9: 14 MMR seems to be the most relevant 
one. Their potential influence on the motion could be investi- 
gated if the initial condition of the system was refined using a 
more accurate model of the dynamics, possibly including rel- 
ativistic effects. However, it would most likely require much 
more precise TOA measurements than currently available. 

These factors allows us to state that the PSR 1257+12 sys- 
tem is orbitally stable over the Gyr time scale. In our exper- 
iments, there are no signs of a potential instability except for 



a very slow divergence of the MEGNO in few of the tests. 
This divergence of MEGNO corresponds to the Lyapunov 
time of about ~ 1 Gyr. We believe that it has only a nu- 
merical character. We are aware that an alternative analytical 
study of the PSR 1257+12 dynamics is possible (also thanks 
to the accurate determination of th e initial condition). Such 
an approa ch has been proposed in Malhotra et al. ( 1982) an d 
Malhotrall ll993al) . Our numerical investigations can certainly 
be treated as a complement to any future analytical studies of 
the system. 

Using the MEGNO analysis, we found dynamical limits on 
the unconstrained elements of the planet A. Due to the desta- 
bilizing effect of the Kozai resonance, the nodal line of this 
planet cannot be separated by more than about ±60° from the 
nodes of B and C. Otherwise, the eccentricity of A would be 
excited to large values permitting close encounters or colli- 
sions with the planet B. It constitutes a strong dynamical ar- 
gument that the orbital plane of the planet A indeed coincides 
with the mean orbital plane of the system. 

Finally, using the MEGNO and the MFT, we investigate the 
dynamics of massless particles in the PSR 1257+12 system in 
the framework of the restricted model. In the numerical ex- 
periments, we find a stable zone between the planets A and 
B extending for initially small eccentricities from 0.19 AU to 
0.25 AU from the pulsar. There are no stable areas between 
the planets B and C. Beyond the orbit of the planet C, the sta- 
ble zone begins already outside of its orbit. We find that the 
massless particles can move on stable orbits under the con- 
dition that their initial eccentricities and semi-major axes are 
located under the collision lines with the planets. The dynam- 
ics of massless particles is basically independent on the their 
initial inclinations. Beyond 1 AU, the motion appears to be 
stable except for the areas of narrow MMRs with the planets 
B and C. It is an encouraging result supporting the search for 
possible small bodies contained in a dust or Kuiper belt type 
disk around the PSR 1257+12 . 
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TABLE 1 

ASTROCENTRIC OSCULATING ORBITAL ELEMENTS OF THE PLANETS IN PSR 1257+12 PLANETARY SYSTEM AT THE 

epoch MJD=49766.50. a 



Planet 


Mass [M e ] 


a [AU] 


e 


i [deg] 


a [deg] 


(0 [deg] 


M [deg] 


A 


0.019 


0.18850 


0.0000 


50.00 


0.00 


0.0 


14.25 


B 


4.250 


0.35952 


0.0186 


53.00 


0.00 


250.4 


5.41 


C 


3.873 


0.46604 


0.0252 


47.00 


3.26 


108.3 


3.66 



a The semi-major axis, a, the eccentricity, e, the inclination, i, the longitude of ascending node, 12, the argument of periastron, 
03, and the mean anomaly, M, of the planets in PS R 1257+12 p lanetary system at the epoch MJD=49766.50 are derived from 
the best-fit orbital parameters by Konacki & Wolszcz aB ETO . The mass of the central star is equal to 1 .4 Mq . 



Pauwels, T. 1983, Celestial Mechanics, 30, 229 Robutel, P. & Laskar, J. 2001, Icarus, 152, 4 

Peale, S. J. 1976, ARA&A, 14, 215 Wisdom, J. & Holman, M. 1991, AJ, 102, 1528 

— . 1993, AJ, 105, 1562 Wolszczan, A. 1994, Science, 264, 538 

Rasio, F. A., Nicholson, P. D., Shapiro, S. L., & Teukolsky, S. A. 1992, Wolszczan, A. & Frail, D. A. 1992, Nature, 355, 145 

Nature, 355, 325 
Rivera, E. J. & Lissauer, J. J. 2001, ApJ, 402, 558 

FIG. 1. — The orbital evolution of the nominal PSR 1257+12 system and its MEGNO signature. Panel (a) is for the semi-major axes. Panel (b) is for the 
orbital inclinations. Panel (c) is for the eccentricities. Panel (d) illustrates the secular apsidal resonance between the planets B and C. Panels (f) and (e) are for the 
time-evolution of MEGNO and its mean value, (Y)(t) (a few representative evolutions for different choices of the initial tangent vector are shown). The orbital 
evolution is given in terms of the heliocentric canonical elements related to the Poincare coordinates. 

FIG. 2. — The dynamics of the PSR 1257+12 system in the 1 Gyr integration. Panel (a) is for the critical argument of the secular apsidal resonance during 
1 Gyr. Panel (b) is for the AMD. Panels (c) and (d) are for the argument 9i = 03a — ®B- Panel (c) is for the first 0.2 Myr while panel (d) is for the end of the 
1 Gyr period. The solid line in these plots denotes the eccentricity of the innermost planet multiplied by 10 4 . Panel (e) shows the SAR between the planets B and 
C in the space of (05^ — i e c)- The same plot for 8j and is shown in panel (f). The orbital evolution is given in terms of the heliocentric canonical elements 
related to the Poincare coordinates. 

FIG. 3. — The dynamical environment of the PSR 1257+12 planets in the space of the semi-major axes. The plots are for the one-dimensional MEGNO scans 
along the semi-major axis of the planet A, B and C. All the other initial elements are fixed at their nominal values (see Table [TV The resolution of the scans is 
5 ■ 10" 6 AU and 10~ 6 AU (the later for the close -up scans for the planets A and C). Labels mark the positions of the MMRs between the planets: the upper plot 
and the second plot for the A-B pair the third, fourth and fifth panel from the top for the B-C pair. The upper plot is a magnification of the scan for the planet A, 
the bottom plot is a magnification of the scan for the planet C. Big dots mark the nominal positions of the planets. 

FIG. 4. — The left column shows the evolution of the critical arguments of the MMRs close to the planet C (see the text for explanation). The panels in the 
right column are for MEGNO. The Lyapunov time about 3000 yr for the C3LB21 MMR is estimated through a linear fit to the MEGNO plot. 

FIG. 5 . — The MEGNO scan along ciq for different masses of the host star. The upper scan is for m, = 0.95Mp Sr and the lower for m* = 1 .05Mp S1 - where Mp Sr 
is the canonical mass of the pulsar, 1 .4 Mq . 

FIG. 6. — The MEGNO stability maps for the configuration given in TablelTl The left panel is for the (ap^.e^-plaiie. The middle panel is for the (ag .eg) -plane. 
The right panel is for the (aQ,eQ) -plane. The position of the nominal system is marked by the two intersecting lines. 

FIG. 7. — The MEGNO stability map for the configuration given in TablelTl The scan is for the (eg -plane. The position of the nominal system is marked 
by the two intersecting lines. 

FIG. 8. — The semi-amplitude 9 max of the SAR in the (<3g,eg ) -plane The position of the nominal system is marked by the two intersecting lines. 

FIG. 9. — The stability map in the (i'Ai^a) -pl ane (the left panel) and the maximal eccentricity of the planet A attained during the integration time (the right 
panel). The resolution of the plot is 1° x 1°. The position of the nominal system is marked by the two intersecting lines. 

FIG. 10. — Kozai resonance for a configuration corresponding to the initial i'a = 50° and Q.^ = —75°. Panel (a) shows the changes of the eccentricity and 
inclination of the planet A relative to the invariant plane of the system. Panel g a is for the argument of periastron, measured with respect to the invariant plane. 
Panel for Y(t) shows the temporal variations of MEGNO and its mean value, (Y). 

FIG. 1 1 . — The stability maps for the massless particles in different regions of the PSR 1257+12 system. The left column is for the region II (between the 
planets A nd B), the middle panel is for the region Ilia (beyond the planet C, up to 1 AU) and the right panel is for the region Illb (up to 3.9 AU). The upper maps 
are for MEGNO, the bottom maps are for the diffusion rate log Co- The resolution of the scans is 0.0008AU x 0.005, 0.0025AU x 0.005 and 0.005AU x 0.1, 
respectively. Collision lines with the planets are marked with white curves. 

FIG. 12. — The MEGNO scan along the semi-major axis ciq of the massless particles moving between the planet A and B (thick line) and their survival times 
over 15 Myr integrations (represented by thin vertical lines). 

FIG. 13. — The stability maps for the massless particles moving in the region Ilia of the PSR 1257+12 system for different initial inclinations. The left column 
is for the the initial inclination in = 75°, the middle column is for to = 87°. The right panels are for eo = and ig 6 [10°, 90°]. The upper maps are for MEGNO, 
the bottom maps are for the diffusion rate logo"o- The resolution of the scans is 0.0023AU x 0.04, 0.0019AU x 0.01, and 0.0019AU x 2°, respectively. 



10 



Gozdziewski, Konacki & Wolszczan 

a) b. a , i3,ic[AlJ| b) i A . i B . i c 

55 .CO 

50 no 



45JQQ 





Q 055 a.75 1 

d?^ - n;,- [deg] 



1 [lfl J jt] 




DQJQ 



:ooi 



:oqg 



025 Q5 a.75 1 



1 [in 3 jt] 



1.5 1[Mjt] 



1.999 I 1 1 1 1 L 

0.5 1 1.5 2 25 1 [Mjt] 



TABLE 2 

Parameters for the (h,k) secular solution obtained by means of the Laplace-Lagrange theory. 

Amplitudes e Pi ; are multiplied by 10 3 . 



Mode [/] 


g [arcsec/yr] 


P [deg] 


«A,i 




e C,i 


Period [yr] 


1 


209.991 


272.28 


-2.075 


21.148 


-20.311 


6171.7 


2 


44.713 


140.17 


-8.146 


0.010 


0.013 


28984.7 


3 


14.061 


333.01 


-6.928 


-7.937 


-7.960 


92169.0 



TABLE 3 

Fundamental frequencies and periods in the PSR 1257+12 system. 



Frequency (period) [i] 1 2 3 



n, [deg/d] 14.249 (25.264 d) 5.410 (66.544 d) 3.665 (98.218 d) 
gi [arcsec/yr] 197.696 (98502.1 yr) 43.881 (29534.0 yr) 13.157 (6555.5 yr) 
.v, [arcsec/yr] -44.023 (29439.3 yr) -203.108 (6380.9 yr) (-) 



t) AMD [Iff^M^lF/d] 
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